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Topological phases are unique states of matter incorporating long-range quan- 
tum entanglement, hosting exotic excitations with fractional quantum statis- 
tics. We report a practical method to identify topological phases in arbitrary 
realistic models by accurately calculating the Topological Entanglement En- 
tropy (TEE) using the Density Matrix Renormalization Group (DMRG). We 
argue that the DMRG algorithm naturally produces a minimally entangled 
state, from amongst the quasi-degenerate ground states in a topological phase. 
This proposal both explains the success of this method, and the absence of 
ground state degeneracy found in prior DMRG sightings of topological phases. 
We demonstrate the effectiveness of the calculational procedure by obtaining 
the TEE for several microscopic models, with an accuracy of order 10~^ when 
the circumference of the cylinder is around ten times the correlation length. 
As an example, we definitively show the ground state of the quantum S =1/2 
antiferromagnet on the kagome lattice is a topological spin liquid, and strongly 
constrain the full identification of this phase of matter. 
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Theory has shown that quantum ground states may exhibit distinct patterns of long-range 
entanglement, which provides the most basic categorization of quantum phases of matter, more 
fundamental than Landau's symmetry breaking paradigm. The simplest and most robust long 
range entangled states, which have a full spectral gap, comprise "topological phases", which 
host topological order. Much recent interest in topological phases is due to the prospect of 
utilizing them to construct an inherently fault tolerant quantum computer (i, 2). Topological 
phases are also of basic scientific interest for their many unique properties, especially their 
ability to support exotic excitations with fractional and even non-abelian quantum statistics. 

Crystalline "Mott" insulators with unpaired electron spins have long been considered likely 
candidates for long range entangled states, epitomized in this context by Anderson's resonating 
valence bond (5) wavefunction for a "quantum spin liquid" (4). Because this particular state is 
non-magnetic, the lack of magnetic order has been widely taken as a definition of a quantum 
spin liquid. However, defining what a quantum spin liquid isn't has little utility, and is especially 
unhelpful in the theoretical search for these phases. Instead, a positive definition of a quantum 
spin liquid which can be tractably tested in realistic models is sorely needed. 

In principle, such a positive definition is provided for topological phases (and hence those 
quantum spin liquids with topological order) by the Topological Entanglement Entropy (TEE) 
introduced by Kitaev-Preskill (5) and Levin- Wen (6). Unfortunately, the formulations in Refs. 
(5, 6) suffer from severe finite-size corrections due to lattice scale effects, greatly hindering 
their application. We report here a practical and extremely simple scheme to numerically cal- 
culate the TEE, and thereby identify topological order. Our method consists simply of using 
the Density Matrix Renormalization Group (7, 8) (DMRG) to calculate the usual entanglement 
entropy for the division of a cylinder into two equal halves by a flat cut, and extracting the TEE 
from its asymptotic, large circumference limit (see below). We argue that this method actually 
works, despite potential complications known theoretically, (9, 10) due to a subtle ground state 
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selection mechanism built into the DMRG algorithm. The approach is tested here on a variety 
of lattice models, and then applied successfully to the physically realistic quantum spin S = ^ 
anti-ferromagnetic Heisenberg J1-J2 model on the kagome lattice. By extracting an accurate 
TEE, we identify a quantum spin liquid state with topological order for the first time in a phys- 
ically realistic S'[/(2)-invariant lattice model. We emphasize that the TEE provides positive, 
"smoking gun" evidence for a topological quantum spin liquid, and excludes any topologically 
trivial states, regardless of possible complex or subtle broken symmetries. The value of the 
TEE also greatly restricts the possible topological quantum field theories which fully describe 
the topological order. We return to this at the end of this paper. 

The TEE is derived from the bipartite von Neumann entanglement entropy, which is defined 
by dividing a system into two subsystems, A and B, which together comprise the full system. 
The entanglement entropy associated to this partition is defined from the reduced density matrix. 
Pa — Tr (|0)(0|), where |0) is a ground state, according to S{A) — — Tr[p^ln(p^)]. It has the 
duality property S{A) = S{B). According to the seminal works of Kitaev-Preskill (5) and 
Levin- Wen (6), the entanglement entropy of a partition of a two dimensional system where A is 
a disk-like region with a smooth boundary (the "entanglement surface") of length i scales as 

S{A) ^ai--f + ---, (1) 

where the ellipsis represents terms that vanish in the limit £ — > 00. The coefficient a, due 
to short distance physics near the boundary, is non-universal. The term 7 is the topological 
entanglement entropy (TEE) — a universal additive constant characterizing the long-range en- 
tanglement in the ground state which can be quantified as 7 = In D, where D is the total 
quantum dimension of the medium (5, 6). 

Note that 7 is sub-dominant to the a£ term, arising from short-range entanglement. As a 
consequence, it is non-trivial to extract. Moreover, for real lattice systems, it is not obvious how 
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to define (. on the lattice, nor is it obvious what qualifies as a "smooth" boundary. These two 
ambiguities are particularly challenging. 

In Refs. (5, (5), complex prescriptions were proposed to remove the short-range contributions 
and extract the TEE from measurements on planar partitions. In our much simpler scheme, we 
study a single partition, defined by a straight cut normal to a cylinder which divides it in half, and 
extract the TEE using Eq. ([T|) with i. = Ly, the circumference of the cylinder. This approach 
minimizes errors due to subtractions of many large numbers, and also minimizes finite size 
corrections due to short-range entanglement, as we now argue. 

For the cylindrical case, we expect such finite size corrections to be of order e'^^l^. In the 
Kitaev-Preskill and Levin- Wen formulations, the corrections are much larger. There, to obtain 
the TEE, the entropy is calculated for several disk-like planar partitions, and comer contribu- 
tions are cancelled by forming a linear combination of the results. However, the complicated 
shape of the planar partitions involved means that the smallest spatial features of the partition 
are several times smaller than the overall system width. For instance, in the Levin- Wen formu- 
lation, the smallest features (size d) are at least four times smaller than the linear width of the 
system assuming periodic boundary conditions, so that L > Ad, a conservative estimate. Cor- 
rections to Eq. ([T]) should be expected to be of order e^'^/^ > e^^^^^^\ Thus to obtain similar 
performance to that of the cylindrical cut, even assuming no additional errors are introduced 
by the subtractions of different entropies, requires a linear system size at least four times larger 
in the Levin-Wen case. This means at least 16 times as many spins, and given the exponential 
growth of the Hilbert space with the number of quantum degrees of freedom, this is a very 
costly increase. Indeed, attempts to implement the Kitaev-Preskill and Levin- Wen protocols in 
simulations have shown them to be very challenging numerically (11, 12). 

A potential complication of our method is that the ground state on a cylinder is expected to 
have a degeneracy in a topological phase in the thermodynamic limit, and the TEE for the cylin- 
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drical cut can depend upon which ground state the TEE is measured in (9, 10). In Ref. (10), it 
has been shown, however, that the TEE for the Z2 spin liquid is bounded above by the universal 
value 7 = In 2, and below by zero. Moreover, in general the universal value is achieved for 
so-called Minimal Entropy States (MES's) (10), which correspond to states in which a quasi- 
particle is definitely contained within the region A (or B). For the Z2 spin liquid, the MES's 
are the states with a definite Z2 magnetic flux through the cyUnder, i.e. the vison or no-vison 
eigenstates. 

We suggest, based on numerical evidence, that the DMRG systematically finds a MES. 
This is perhaps natural since the DMRG prefers low entanglement states (8). Note though 
that the absolute ground state of a finite system is dependent upon microscopic details, and is 
expected to vary with the aspect ratio (L^/Ly) of the cyUnder (13). For a physical Hamiltonian 
without fine tuning, the absolute ground state becomes a MES in the "long" cylinder limit, 
where L^/ Ly is larger than some critical value (which depends on microscopic details, but is 
order one generically) (13). However, we contend that the DMRG preferentially finds the MES 
even when it is not the absolute ground state. Evidence for this is given below in the toric code 
model, where the MES can be explicitly identified. The fact that we obtain the universal value 
of the TEE, independent of the system's aspect ratio, for several other models, also supports 
this conclusion. 

We turn now to the toric code model, which is well known - see the Supplementary material 
for details of the definition. It can be considered as a model of fluctuating discrete "electric" 
and "magnetic" fields. To observe the ground state selection, we first consider an applied field 
h — hx which is purely electric, — 0. Then the operator G (defined in the Supp. Mat.), 
which measures the parity of the number of electric field loops winding around the cylinder, 
commutes with the Hamiltonian, so the energy eigenstates must also be eigenstates of G = ±1. 
Topological order implies that there are two such states with G — +1 and G — —1, with 
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exponentially close energies. The MES's, however, are not G (or energy) eigenstates, but rather 
the superpositions |±) = (|G = I) ±\G = -l))/v^, for which (±|G|±) = 0. The |±) 
states correspond to states with or without magnetic flux through the cylinder. Measurements 
of (G) and 5*, Figure [T| show that the DMRG preferentially selects a MES for larger systems, 
and that the number of states m necessary to converge to the absolute ground state (with larger 
entanglement and zero TEE - see Figure |2j5) grows very rapidly with system size. 

The origin of the topological contributions to the entanglement entropy sheds light on this 
behavior in the Z2 case of interest. First, there is a reduction of entropy due to the constraint 
that electric field loops always cross the entanglement surface an even number of times. This 
reduction is precisely the TEE, and this physics is included once entanglement on the scale of Ly 
is taken into account. Second, in the case where the absolute ground state is not a MES, there is 
an increase of entropy due to the global constraint on the number of electric field lines winding 
the entire cylinder. To take this into account, the DMRG must fully converge the entanglement 
of the opposite ends of the system, which are extremely far separated on the "snaking" DMRG 
path. This global entanglement does not converge for larger systems, in which case the DMRG 
produces states described by a Schmidt decomposition in which the left and right halves of the 
system have uncorrelated electric field winding parities. Such a state is a MES. 

We next consider the toric code model in symmetrically applied fields, = = h, for 
which the absolute ground state is not obvious. Figure [2] shows the entanglement entropy in this 
case. This model was previously shown (14) to have a quantum phase transition between the 
Z2 phase for /i < /ic 0.34 and a trivial phase for h > h^. The extrapolated TEE following our 
protocol indeed very well approximates the universal value 7 = ln2 = 0.69314 ... for /i < he, 
even for h = 0.3, relatively close to the quantum phase transition, we obtain 7 = 0.691(4), 
which is accurate to a fraction of a percent. For h > we obtain 7 = 0, as expected, with a 
numerical uncertainty of order 10^^. Similar results are obtained for a variety of aspect ratios 
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Figure 1: Evidence that the DMRG favors MES's. In (a) the electric field parity {G) and (b) the 
entanglement entropy are plotted versus number of states m for the toric code with = 0.3 
and hz = 0, for several system sizes. We see that for fixed small system size, at smaller m the 
average parity (G) is approximately zero and the entanglement is reduced, while at large m a 
definite parity eigenstate is found with G = 1, and the entanglement is increased. The jumps 
in the two quantities coincide, signaling a transition from a MES to an absolute Hamiltonian 
eigenstate. The number of states m needed to converge to the absolute ground state increases 
rapidly with L^. For larger systems than shown, a MES with (G) ^ Ois found for all accessible 
values of m. 
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and values of the perturbing fields. 

We apply the method to the spin- 1/2 Heisenberg model on the kagome lattice, for which 
compelling but indirect evidence for a gapped quantum spin liquid has been recently obtained by 
extensive DMRG studies (15-17). We consider the model with both first and second-neighbor 
interactions, whose Hamiltonian is 



where Sj is the spin operator on site i, and {ij) ({{ij))) denotes the nearest neighbors (next 
nearest neighbors). In the numerical simulation, we set Ji = 1 as the unit of energy. The most 
recent DMRG studies (17) show that the J2 = point is near the edge of a substantial spin 
liquid phase centered near J2 = 0.05 — 0.15. 

We take the kagome lattice with periodic boundary conditions along a bond direction, drawn 
vertically in the inset of Figure [3j and the unit of length equal to the nearest-neighbor distance. 
The results for the entanglement entropy for J2 = 0.10 and 0.15 are shown in Figure |3] with 
correlation length around one-lattice spacing for both spin- spin and dimer-dimer correlation 
functions. We see that a linear fit using data for = 4 ~ 12 using Eq.Q gives 7 = 0.698(8) 
at J2 = 0.10 and 7 = 0.694(6) at J2 = 0.15, both within one percent of In 2 = 0.693. This 
proves definitively that this phase is a topological spin liquid, and determines the quantum 
dimension D = 2, very consistent with a Z2 state. 

We have shown that the TEE can be calculated to an accuracy of order 10^'^ when Ly is 
~10 times the correlation length (see the Supporting Online Material for some additional tests). 
Our result provides a "smoking gun" test for a topological spin liquid. It also explains the 
puzzling absence of topological degeneracy in recent DMRG results which otherwise support 
a Z2 spin liquid state (13, 16), since we have shown that the DMRG is systematically biased to 
find just one of the ground states. The TEE does not fully determine the nature of the topological 
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Figure 2: The von Neumann entropy S{Ly) for the toric code model in magnetic fields. In (a), 
we show S{Ly) with = 4 ~ 16 at = oo for symmetric magnetic fields at = = h = 
0.2, 0.3 and 0.4. By fitting S{Ly) = aLy - 7, we get 7 = 0.693(1), 0.691(4) and 0.001(5), 
respectively. In (b), we consider the pure electric case, = 0.3, = 0, and compare S{Ly) in 
the MES obtained in the large limit (black squares) to that of the absolute ground state from 
systems of dimensions Lj. x Ly = 20 x A, 2A x 6, 24 x 8, 24 x 10 (red circles). Extrapolation 
shows that the MES has the universal TEE, while the absolute ground state has zero TEE. 
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Figure 3: The entanglement entropy S{Ly) of the kagome J1-J2 model in Eq.([2]), with Ly = 
4 ~ 12 at L.^ = 00. By fitting S{Ly) = aLy - 7, we get 7 = 0.698(8) at J2 = 0.10, and 
7 = 0.694(6) at J2 = 0.15. Inset: kagome lattice with = 12 and Ly = 8. 

phase. Fortunately, for a given D, there are only finitely many distinct topological phases, and 
for small values of D, a complete classification of all topological phases is known (18). Other 
constraints such as time -reversal symmetry (if present) further constrain the possible topological 
order. For example, there are only two time-reversal invariant phases consistent with D = 2, 
found here for the kagome Heisenberg model: the Z2 phase, and a doubled semion phase. 
It will be interesting to develop methods to distinguish these in the future, and to calculate 
the topological ground state splitting. Identifying topological order by combining theoretical 
classification results with numerical simulation is a major step in the development of a post- 
Landau paradigm for classifying quantum phases of matter. 
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Supporting Online Material 

Materials and Methods 

Here we test our method on a variety of lattice models whose topological order is known. 

A Toric-code model in magnetic fields 

The toric code model (i) with an applied magnetic field is given by 

i7 = -J,5^A,- Jp5^5p-/i,5^af-/i,5^af, (SI) 

s p i i 

where and are Pauli matrices, and Ag = HjgsO'f , Bp = Ui^pa-. Subscripts s and p refer 
to, respectively, vertices and plaquettes of a square lattice, whereas i runs over all bonds where 
spin degrees of freedom are located. Without magnetic field, i.e., = = 0, the pure toric 
code model can be solved exactly (i), and the ground state has Z2 topological order with total 
quantum dimension D = 2. On the torus the ground state is 4-fold degenerate. All elementary 
excitations are gapped and characterized by eigenvalues As = — 1 (a Z2 charge on site s) and 
Bp = —1 (a. Z2 vortex on plaquette p). When turning on the magnetic field, the model cannot 
be solved exactly anymore. Previous studies (14, 19, 20) show that the Z2 topological phase 
remains stable and robust until the magnetic fields are large enough, where the system crosses 
the transition from the topological phase to the trivial one. Specifically, such a phase transition 
takes place at the critical magnetic field h^, = 0.34 along the symmetric line = = h. 

For the DMRG simulation, we consider an equivalent square lattice, where the spin oper- 
ators cr^ and cr^ sit on the sites instead of the bonds. Therefore, the star operator Ag and the 
plaquette operator Bp of the original lattice now sit on alternating plaquettes in the equivalent 
square lattice, as shown in Figure [ST] labeled as S and P, respectively. Note that on this equiv- 
alent square lattice, there are an even number of dangling spins within each plaquette at the 
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Figure SI: The von Neumann entropy for the pure toric code model and that in magnetic fields 
in Eq.(Sl I. The von Neumann entropy S{Ly) for the pure Toric-Code model with Lj^ = 4 ~ 20 
at Lx = oo. By fitting S{Ly) = aLy — 7, we get 7 = 0.693147(1). Inset: (Upper) the 
entanglement entropy S{Lx) as a function of for different L^; and (Lower) Square lattice 
with = 10 and Ly = 6. Here S represents the star operator Ag, while P represents the 
plaquette operator Bp. 



open edges. For the pure toric-code model with cylinder boundary condition, the first 2^y^'^~^ 
eigenvalues of the reduced density matrix are degenerate and equal to 1/2^^'/^^^, while all 
the other eigenvalues are zero. This allows us to study a quite large system with width up to 
Ly = 20 easily. As shown in Figure SI after fitting the entanglement entropy using Eq.Q for 
Lj, = 00, we get a nonzero topological entanglement entropy 7 = 0.693147(1), which is equal 
to the expected value 7 = In 2 with amazingly small numerical error 10^^. 

After turning on the magnetic field, the degeneracy of entanglement spectrum is lifted, and 
the correlation length ^ becomes finite. Our results show that even very close to the phase 
transition point (e.g., he = 0.34 along the symmetric line), we can still get a very accurate TEE 
7. For example, the correlation length ^ ~ 1 lattice spacing ath = 0.30, the resulted topological 
entanglement entropy 7 = 0.691(4) is still quite accurate with an error around 10^^. These 
results show that a nonzero TEE is obtained throughout the topologically ordered phase. On 
the contrary, the topological entanglement entropy 7 is zero in the trivial phase where h > he. 
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For example, 7 = 0.001(5) at h — 0.40. Therefore, our method allows us to unambiguously 
extract the non-zero topological entanglement entropy 7 if and only if the toric-code model is 
in a topologically ordered phase. 

In the case of a purely electric perturbation, hz — 0, hx — h 0, two loop operators 
commute with H. Specifically, these are 

G = G, = n<., (S2) 

x=l 

Gx = (S3) 

y=l 

In the low energy sector where Ag — +1 for all s, Gy is independent of y and Gx is independent 
of X. By construction, G and Gx have eigenvalues ±1. The operator Gx probes the presence or 
absence of an electric particle at the end of the cylinder. This degree of freedom is not associated 
to the ground state degeneracy, and indeed we find Gx — +1 always in our numerics. The 
operator G counts the parity of the number of electric flux lines winding around the cylinder, and 
does operate in the topologically degenerate subspace. Physically, eigenstates of G are equal 
weight superpositions of the vison and no-vison eigenstates, which are the MES, as discussed 
in the main text. For hz — 0, the energy eigenstates must also be eigenstates of G, and the 
splitting between them is expected to be expontially small in Lx. 

B Fractional quantum Hall model 

We next consider the so-called Haldane model (21) on the honeycomb lattice filled with hard- 
core bosons: 

H ^ -t' J2 [^'^r-e^^^'^ + H.c] (S4) 

{{rr'» 

- t J2 [bl'br + H.c] - t" ^ ^blK + H.c] , 



{rr') (((rr'»> 
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where 6j creates a hard-core boson at site r, (■ ■ ■), ((■ ■ ■)), and (((■ ■ ■))) denote the nearest- 
neighbor, the next-nearest-neighbor, and the next-next-nearest-neighbor pairs of sites, respec- 
tively. In Ref. (22), the authors have systematically studied this model using exact diagonaliza- 
tion, providing convincing evidence showing that the ground state (with parameters t' = 0.6t, 
t" = — 0.58t, and = OAtt) is a 1/2 bosonic FQH state with two-fold ground state degeneracy 
on the torus. Such a 1/2 FQH state has nontrivial semion topological order (5, 6), with total 
quantum dimension D = \pl. 

For the numerical simulation, we consider a honeycomb lattice with length vectors Liai 



and 1^232 as shown in the inset (a) of Figure S2 Here ai = (v^, 0) and 02 = |) are 



two primitive vectors of the unit cell which includes two sites of the lattice. The total number 
of sites is = 2 X Li X L2, with Li x L2 unit cells. Note that the corresponding system 
width Ly = 2L2, and system length = 2Li. Unambiguously, extrapolation from the data 
for Ly < 20 using Eq.Q shows that we can get a nonzero constant topological entanglement 
entropy 7 = 0.349(5). This is equal to \n(\/2) = 0.347 within the numerical error, showing 
that our method can also be used to study chiral topological states as well. 

C Transverse-field Ising model 

The models studied above have topologically ordered ground states, for which our method in- 
deed gives us non-zero topological entanglement entropy with high accuracy. Now, we show 
that for a topologically trivial phase, our method unambiguously gives zero topological entan- 
glement entropy as well. To show this, consider the well-known transverse field quantum Ising 
model 

H = -J2^t^j-hJ2^!, (S5) 
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Figure S2: The entanglement entropy S{Ly) of the Honeycomb Haldane model in Eq.(S5), 
with = 6 ~ 20 at = oo. By fitting S{Ly) = aLy — 7, we get 7 = 0.349(5). Inset: (a) 
Honeycomb lattice with Li = 6 and L2 = 4. Here the system width Ly = 2L2, and system 
length Lx = 2Li. (b) The entanglement entropy S{Lx) as a function of Lj. for different Ly. 



where and a- are Pauli matrices on site i. This model is known to have a topological trivial 
ground state in all magnetic fields, and a second order phase transition at critical field he = 
3.044 (23). As shown in Figure [S3]^a), our method produces very accurate results showing 
7 = even very close to the phase transition point. For example, at h = 3.1, we obtain 
7 = 0.0014(5), which is zero within the numerical error, despite the longish correlation length 
^ ~ 4 for < S'^S^ > and ^ ~ 1 for < S^'S'-' > in this case. 

D Coupled spin-dimer model 



Another well-known model with a topologically trivial ground state is the coupled spin-dimer 
model. 



H = ^Si-Sj + g ^ Si - Sj, 

{ij)eA {ij)&B 



(S6) 



where Sj is the spin-^ operator on site i on the square lattice shown in the inset of Figure S3 'b) 



with A hnks forming decoupled dimers while B links couple the dimers. The ground state of 
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Figure S3: The entanglement entropy of the transverse field quantum Ising model and coupled 
spin-dimer model, (a) The entanglement entropy S{Ly) of the transverse field quantum Ising 
model in Eq.(S5) with Ly = 4 ~ 24 at L^; = oo, in different magnetic field h = 3.1, 3.5 and 
4.0. By fitting S{Ly) = aLy - 7, we get 7 = 0.0014(5), 0.0004(4) and 0.0001(2), respectively. 
Inset: The entanglement entropy S{Lx) as a function of for different Ly at h = 4.0. (b) The 
entanglement entropy S{Ly) of the coupled spin-dimer model in Eq.(S6) with Ly = A 16 
at L3. = 00, at different coupling g = 0.5, 0.4 and 0.3. By fitting S{Ly) = aLy — 7, we get 
7 = 0.006(6), 0.002(1) and 0.0008(9), respectively. Inset: The coupled spin-dimer model, with 
spin (5* = ^) on the sites, the A links are shown as full lines, and the B links as dashed lines. 
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Eq.(S6) depends only on the dimensionless coupling g. It is known that there is an gapped 
dimerized phase for g < Qc = 0.52 {24, 25), at which point a second order phase transition oc- 



curs. Unlike in the transverse field quantum Ising model in Eq.(S5), spin rotational symmetry 
is preserved in this model, although the lattice translational symmetry is explicitly broken. Un- 



ambiguously, as shown in Figure S3 ^b), our method once more produces very accurate results 
showing 7 = even quite close to the phase transition point. For example, 7 = 0.006(6) at 
g = 0.50, which is zero within the numerical error. 
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